Integrating Primate Data into Analysis

Setup

Load libraries

library(ggplot2)
library(tidyr)
library(dplyr)
library(Matrix)
library(Seurat)
library(cowplot)
library(patchwork)

# parallelization
library(future)
options(future.globals.maxSize= +Inf)
plan()

sequential:
- args: function (expr, envir = parent.frame(), substitute = TRUE, lazy = FALSE, seed = NULL, globals = TRUE, local = TRUE, earlySignal = FALSE, label = NULL, ...)
- tweaked: FALSE
- call: NULL

Process Human Data

import_remote_data <- function(file_url, type = "table", header = FALSE) {
  con <- gzcon(url(file_url))
  txt <- readLines(con)
  if (type == "MM") { return (readMM(textConnection(txt))) }
  if (type == "table") { return (read.table(textConnection(txt), header = header)) }
}
count_matrix_URL <- "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE137nnn/GSE137537/suppl/GSE137537_counts.mtx.gz"
gene_names_URL <- "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE137nnn/GSE137537/suppl/GSE137537_gene_names.txt.gz"
sample_annotations_URL <- "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE137nnn/GSE137537/suppl/GSE137537_sample_annotations.tsv.gz"

human.count_matrix <- as.matrix(import_remote_data(count_matrix_URL, type = "MM"))
human.gene_names <- import_remote_data(gene_names_URL, type = "table")
human.sample_annotations <- import_remote_data(sample_annotations_URL, type = "table", header = TRUE)

human_ret_seurat

An object of class Seurat 
19712 features across 20091 samples within 1 assay 
Active assay: RNA (19712 features, 0 variable features)

Process Mouse Data

mouse.data <- Read10X(data.dir = "filtered_feature_bc_matrix")
dimnames(mouse.data)[[1]] <- tolower(dimnames(mouse.data)[[1]])
dimnames(mouse.data)[[2]] <- tolower(dimnames(mouse.data)[[2]])
mouse_ret_seurat <- CreateSeuratObject(counts = mouse.data, 
                                       project = "mouse_ret", 
                                       min.cells = 3, 
                                       min.features = 200)
mouse_ret_seurat

An object of class Seurat 
16424 features across 4510 samples within 1 assay 
Active assay: RNA (16424 features, 0 variable features)

Process Primate Data

url=https://ftp.ncbi.nlm.nih.gov/geo/series/GSE118nnn/GSE118546/suppl/GSE118546_macaque_fovea_all_10X_Jan2018.Rdata.gz
wget $url -O primate_data/GSE118546_macaque_fovea_all_10X_Jan2018.Rdata.gz
gunzip primate_data/*

install.packages( c('devtools', 'roxygen2') )
library(devtools)
library(roxygen2)
install_github( 'hb-gitified/cellrangerRkit',
                auth_token = 'your_token' )

macaque_fovea_seurat

An object of class Seurat 
30039 features across 111993 samples within 1 assay 
Active assay: RNA (30039 features, 0 variable features)

Cleanup

rm(human.count_matrix, human.gene_names, human.sample_annotations)
rm(count_matrix_URL, gene_names_URL, sample_annotations_URL, import_remote_data)
rm(mouse.data)
rm(Count.mat_fovea, macaque_fovea)

Combine

# combine
ret.list <- list(human = human_ret_seurat, mouse = mouse_ret_seurat, macaque = macaque_fovea_seurat)

# preprocess
ret.list <- lapply(X = ret.list, FUN = function(x) {
    x <- NormalizeData(x, verbose = FALSE)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000, verbose = FALSE)
})

# cleanup
rm(human_ret_seurat, mouse_ret_seurat, macaque_fovea_seurat)

Integration

plan("multiprocess", workers = 4)
ret.anchors <- FindIntegrationAnchors(object.list = ret.list, dims = 1:50,  anchor.features = 1000)
plan("multiprocess", workers = 1)
ret.combined <- IntegrateData(anchorset = ret.anchors, dims = 1:50)

Integrated Analysis

plan("multiprocess", workers = 4)

DefaultAssay(ret.combined) <- "integrated"

# Run the standard workflow for visualization and clustering
ret.combined <- ScaleData(ret.combined, verbose = FALSE)
ret.combined <- RunPCA(ret.combined, npcs = 50, verbose = FALSE)
# t-SNE and Clustering
ret.combined <- RunUMAP(ret.combined, reduction = "pca", dims = 1:35)
ret.combined <- FindNeighbors(ret.combined, reduction = "pca", dims = 1:35)
ret.combined <- FindClusters(ret.combined, resolution = 0.075)

UMAP Visualization

DimPlot(ret.combined, reduction = "umap", group.by = "orig.ident")

DimPlot(ret.combined, reduction = "umap", label = TRUE)

DimPlot(ret.combined, reduction = "umap", split.by = "orig.ident", ncol = 1)

Identify Clusters with Canonical Markers

DefaultAssay(ret.combined) <- "RNA"

features <- tolower(c("Pde6a","Gnat2","Nefl","Camk2b","Thy1","Gad1","Slc6a9",
                      "Pcsk6","Trpm1","Sept4","Glul","Arr3","C1qa","Tm4sf1", "Mgp"))

FeaturePlot(object = ret.combined, 
            features = features, 
            pt.size = 0.1,
            cols = c("lightgrey", "#F26969"),
            min.cutoff = "q9",
            combine = TRUE) & NoLegend() & NoAxes()

# for(i in 1:length(p)) {
#   p[[i]] <- p[[i]] + NoLegend() + NoAxes()
# }
# 
# cowplot::plot_grid(plotlist = p, ncol=3)

• Rod : pde6a
• AC (amacrine cell) : gad1, slc6a9
• MG (Müller glia) : glul
• BC (bipolar cell) : Trpm, camk2b
• CC (cone cell) : gnat2, arr3
• RGC (retinal ganglial cell) : nefl, thy1
• VC (vascular cell) : mgp, tm4sf1
• M (microglia) : c1qa
• HC (horizontal cell) : sept4

Markers were determined from this paper and other sources.

ret.combined <- RenameIdents(ret.combined, `0` = "MG", `1` = "Rod", `2` = "RGC", 
    `3` = "RGC", `4` = "BC", `5` = "CC", `6` = "BC", `7` = "AC", `8` = "BC", `9` = "RGC", 
    `10` = "RGC", `11`= "HC", `12` = "MG", `13` = "VC", `14` = "RGC", `15` = "RGC", `16` = "M", `17` = "RGC")

DimPlot(ret.combined, label = TRUE)

Find Differentially Expressed Genes

ret.combined$celltype.organism <- paste(Idents(ret.combined), ret.combined$orig.ident, sep = "_")
ret.combined$celltype <- Idents(ret.combined)
Idents(ret.combined) <- "celltype.organism"

cells.diffgenes <- as.list(cells.types)
cells.diffgenes <- lapply(cells.diffgenes, FUN = function(x) {
  lab_human <- sprintf("%s_human_ret", x)
  lab_mouse <- sprintf("%s_mouse_ret", x)
  return(FindMarkers(ret.combined, ident.1 = lab_human, ident.2 = lab_mouse, verbose = FALSE))
})

Tables with the most differentially expressed genes in each cell subtype:

for(i in seq_along(cells.diffgenes)) {
  print(knitr::kable(head(cells.diffgenes[[i]]),caption=cells.types[[i]]))
}

Rod

	p_val	avg_logFC	pct.1	pct.2	p_val_adj
ckb	0	1.4493770	0.918	0.724	0
hsp90aa1	0	1.3457646	0.854	0.627	0
nrl	0	1.3140138	0.874	0.635	0
0610009b22rik	0	-0.6622860	0.000	0.130	0
gm17018	0	-0.6831275	0.000	0.130	0
spata1	0	-0.6929677	0.000	0.132	0

BC

	p_val	avg_logFC	pct.1	pct.2	p_val_adj
neat1	0	3.086391	0.793	0.064	0
mtch1	0	-1.305054	0.000	0.459	0
selenom	0	-1.338108	0.000	0.480	0
araf	0	-1.342891	0.013	0.494	0
klc3	0	-1.424615	0.002	0.500	0
pea15a	0	-1.427543	0.000	0.500	0

MG

	p_val	avg_logFC	pct.1	pct.2	p_val_adj
tf	0	5.089073	0.962	0.000	0
spp1	0	3.879036	0.847	0.003	0
crabp1	0	3.865908	0.876	0.028	0
gpx3	0	3.736219	0.869	0.052	0
ftl	0	3.672007	0.877	0.000	0
actg1	0	3.639157	0.905	0.026	0

RGC

	p_val	avg_logFC	pct.1	pct.2	p_val_adj
mt-nd4	0	-5.434720	0	1	2e-06
mt-nd5	0	-4.555808	0	1	2e-06
mt-co1	0	-4.634061	0	1	2e-06
malat1	0	-5.358199	0	1	2e-06
mt-nd1	0	-5.600755	0	1	2e-06
mt-nd2	0	-5.700498	0	1	2e-06

CC

	p_val	avg_logFC	pct.1	pct.2	p_val_adj
gm42418	0	-5.444663	0	1	0
malat1	0	-5.893437	0	1	0
mt-cytb	0	-6.148057	0	1	0
mt-co1	0	-4.052888	0	1	0
mt-nd5	0	-4.170730	0	1	0
mt-nd1	0	-4.734329	0	1	0

AC

	p_val	avg_logFC	pct.1	pct.2	p_val_adj
mt-nd5	0	-3.999182	0	1	0
gm42418	0	-5.868449	0	1	0
mt-co1	0	-4.023145	0	1	0
mt-nd4	0	-5.035458	0	1	0
mt-nd1	0	-5.173153	0	1	0
mt-nd2	0	-5.370662	0	1	0

VC

	p_val	avg_logFC	pct.1	pct.2	p_val_adj
hla-b	0	3.429125	0.884	0.00	0
rps3a	0	2.938618	0.826	0.00	0
hla-e	0	3.220179	0.826	0.00	0
hla-a	0	3.030967	0.812	0.00	0
hla-c	0	2.913989	0.797	0.00	0
a2m	0	3.322554	0.797	0.01	0

HC

	p_val	avg_logFC	pct.1	pct.2	p_val_adj
mt-nd5	0	-4.723687	0	1	0
mt-co1	0	-4.915798	0	1	0
mt-nd4	0	-5.757623	0	1	0
mt-nd1	0	-5.944134	0	1	0
gm42418	0	-6.046691	0	1	0
mt-nd2	0	-6.094723	0	1	0

M

	p_val	avg_logFC	pct.1	pct.2	p_val_adj
ftl	0	4.612295	0.98	0	0
hla-dra	0	4.664191	0.94	0	0
hla-a	0	2.899218	0.94	0	0
hla-drb1	0	4.096260	0.92	0	0
rps3a	0	3.397725	0.92	0	0
hla-b	0	3.163581	0.92	0	0

Save as csv files

for(i in seq_along(cells.diffgenes)) {
  write.csv(cells.diffgenes[[i]], sprintf("results/%d_%s.csv", i, cells.types[[i]]))
}

genes_to_plot <- 3
for (i in seq_along(cells.types)) {
  print(FeaturePlot(object = ret.combined, 
              features = rownames(cells.diffgenes[[i]])[1:genes_to_plot], 
              split.by = "orig.ident", 
              max.cutoff = 3, 
              cols = c("grey", "red"),
              pt.size = 0.07,
              combine = TRUE,
              label.size = 0.5
              ) + plot_annotation(title = cells.types[[i]]) & NoLegend() & NoAxes()
        )
}

Check cell proportion for each species:

knitr::kable(prop.table(x = table(Idents(ret.combined), ret.combined@meta.data$orig.ident), margin = 2))

	human_ret	macaque_fovea	mouse_ret
0	0.2627047	0.1535810	0.2875831
1	0.5498980	0.0758172	0.3164080
2	0.0001493	0.1855205	0.0002217
3	0.0009955	0.1458216	0.0035477
4	0.0531581	0.0808176	0.0838137
5	0.0114977	0.0783888	0.0388027
6	0.0406650	0.0552267	0.0569845
7	0.0187148	0.0577625	0.0343681
8	0.0484794	0.0443242	0.0917960
9	0.0001493	0.0342075	0.0000000
10	0.0000498	0.0232247	0.0004435
11	0.0076154	0.0156081	0.0152993
12	0.0000000	0.0117775	0.0000000
13	0.0034344	0.0089827	0.0436807
14	0.0000000	0.0109203	0.0000000
15	0.0000000	0.0101435	0.0008869
16	0.0024887	0.0039645	0.0261641
17	0.0000000	0.0039110	0.0000000